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Abstract 

In  many  engineering  applications,  including  surveillance,  guidance,  or  navigation,  single  stand-alone  sen¬ 
sors  or  sensor  networks  are  used  for  collecting  information  on  time  varying  quantities  of  interest,  such  as 
kinematical  characteristics  and  measured  attributes  of  moving  or  stationary  objects  of  interest  (e.g.  maneu¬ 
vering  air  targets,  ground  moving  vehicles,  or  stationary  movers  such  as  a  rotating  antennas). 

More  strictly  speaking,  in  these  applications  the  state  vectors  of  stochastically  moving  objects  are  to  be 
estimated  from  a  series  of  sensor  data  sets,  also  called  scans  or  data  frames.  The  individual  measurements 
are  produced  by  the  sensors  at  discrete  instants  of  time,  being  referred  to  as  scan  or  frame  time,  target  revisit 
time,  or  data  innovation  time.  These  output  data  (sensor  reports,  observations,  returns,  hits,  plots)  typically 
result  from  complex  estimation  procedures  themselves  characterizing  particula  r  waveform  parameters  of  the 
received  sensor  signals  ( signal  processing). 

In  case  of  moving  point-source  objects  or  small  extended  objects,  i.e.  typical  radar  targets,  often  rela¬ 
tively  simple  statistical  models  can  be  derived  from  basic  physical  laws  describing  their  temporal  behavior 
and  thus  defining  the  underlying  dynamical  system.  In  addition,  appropriate  sensor  models  are  available 
or  can  be  constructed,  which  characterize  the  statistical  properties  of  the  produced  sensor  data  sufficiently 
correct. 

As  an  introduction  to  target  tracking  and  data  fusion  applications  characteristic  problems  occurring  in 
typical  radar  applications  are  presented;  key  ideas  relevant  for  their  solution  are  discussed. 
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1  Adaptive  Tracking  for  Phased-array  Radar 

Modern  phased-array  radars  call  for  tracking  algorithms  that  efficiently  exploit  all  degrees  of  freedom  avail¬ 
able,  which  arc  variable  over  a  wide  range  and  can  be  chosen  individually  for  each  track.  This  is  particularly 
true  for  military  air  situations  where  both,  agile  targets  and  targets  significantly  differing  in  their  radar  cross 
section,  must  be  taken  into  account.  Unless  properly  handled,  such  situations  can  be  highly  allocation  time 
and  energy  consuming.  In  this  context  adaptive  techniques  for  combined  tracking  and  sensor  management 
are  discussed,  i.e.  control  of  data  innovation  intervals,  radar  beam  positioning,  and  transmitted  energy.  By 
efficient  exploitation  of  the  limited  resources  the  total  surveillance  performance  of  the  sensor  system  can  be 
improved. 


1.1  General  Overview 

In  military  applications  often  distinct  maneuvering  phases  exist  as  even  agile  targets  do  not  permanently 
maneuver.  Nevertheless,  abrupt  transitions  to  high-g  turns  can  well  occur.  Allocation  time  and  energy 
savings  are  thus  to  be  expected  if  adaptive  target  dynamics  models  are  used[2,  4,  1 1,  20],  such  as  Interacting 
Multiple  Model  (IMM).  Besides  their  kinematical  characteristics,  the  mean  radar  cross  section  (RCS)  of  the 
targets  is  usually  unknown  and  variable  over  a  wide  range.  By  processing  of  signal  amplitude  information, 
however,  the  energy  spent  for  track  maintenance  can  be  adapted  to  the  actual  target  strength.  Also  by  this 
measure  the  total  sensor  load  can  significantly  be  reduced. 

Due  to  the  locally  confined  target  illumination  by  the  pencil-beam  of  a  phased-array,  abrupt  transitions 
into  maneuvering  flight  phases  are  critical  as,  in  contrast  to  conventional  radar,  a  periodic  target  illumination 
is  no  longer  guaranteed.  Any  track  reiinitiation  is  thus  highly  allocation  time  and  energy  consuming  and 
in  addition  locks  the  sensor  for  other  tasks  (e.g.  weapon  guidance).  This  calls  for  intelligent  algorithms  for 
beam  positioning  and  local  search  [8]  that  are  crucial  for  phased-array  tracking. 

For  radars  scanning  the  field  of  view  at  a  constant  rate  (TWS:  Track- while-scan),  the  BAYESian  track¬ 
ing  techniques,  discussed  in  the  first  talk  provide  an  iterative  updating  scheme  for  conditional  probability 
densities  of  the  target  state  given  all  sensor  data  and  a  priori  information  available.  In  those  applications 
data  acquisition  and  tracking  are  completely  decoupled.  For  phased-array  radar,  however,  the  current  signal- 
to-noise  ratio  of  the  target  (i.e.  the  detection  probability)  strongly  depends  on  the  correct  positioning  of  the 
pencil-beam  which  is  now  taken  into  the  responsibility  of  the  tracking  system.  Thus  sensor  control  and 
data  processing  are  closely  interrelated.  This  basically  local  character  of  the  tracking  process  constitutes  the 
principal  difference  between  phased-array  and  TWS  applications  from  a  tracking  point  of  view.  By  using  a 
suitable  sensor  model,  however,  this  fact  can  be  incorporated  into  the  BAYESian  formalism. 

Figure  1  provides  a  simplified  schematic  overview  of  phased-array  radar  function  control  with  the  corre¬ 
sponding  flow  of  information. 


1.2  Modeling  Assumptions 

Resource  management  for  a  multi-functional  phased-array  radar  certainly  depends  on  the  particular  appli¬ 
cation  considered.  We  here  discuss  track  maintenance  for  ground-based  air  surveillance  while  minimizing 
the  allocation  time  and  energy  required.  The  track  accuracy  is  important  only  in  so  far  as  stable  tracks  arc 
guaranteed.  Track  initiation  or  implementation  issues  are  not  addressed  here.  To  make  the  benefits  of  IMM 
modeling  and  amplitude  information  clearly  visible,  false  detections  (clutter,  ECM),  data  association  con¬ 
flicts,  or  possibly  unresolved  measurements  were  excluded  from  the  discussion.  Nevertheless,  their  impact 
might  well  be  incorporated  into  the  general  BAYESian  framework  sketched  in  the  first  talk  (see  also  [7]). 
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Figure  1 :  Phased- Array  Radar  Function  Control  (Flow  of  Information) 

1.2.1  Resource  Allocation 

In  phased-array  tracking  additional  sensor  information  can  be  acquired  when  needed.  Before  each  ‘radar 
resource  allocation’  [4]  certain  radar  parameters  must  be  selected  by  the  tracking  system  depending  on  the 
current  lack  of  information.  We  here  consider  the  target  revisit  time  tk,  the  current  beam  position  and 
the  transmitted  energy  per  dwell  e*.  Other  parameters,  such  as  detection  threshold  Ad  or  beam  width  B, 
are  assumed  to  be  constant.  After  processing  the  skin  echo,  the  resource  allocation  R k  at  time  tk  results  in 
measurements  of  direction  cosines  and  target  range,  z k  =  (uk ,  W  AJ,  along  with  the  signal  amplitude  a*. 
Possibly  a  single  dwell  is  not  sufficient  for  target  detection  and  a  subsequent  fine  localization.  Let  n/,k  denote 
the  number  of  dwells  for  a  successful  detection  along  with  the  corresponding  beam  positions  Bk  =  { bj,. } .  ** . 
Each  allocation  is  thus  characterized  by  the  tuple  Rk  =  (tk,  Bk,nbk,  ek,Zk,  ak )•  The  sequence  of  successive 
allocations  is  denoted  by  Zk  =  {Rk,  Zk~l }. 

1.2.2  RCS  Fluctuations 

The  instantaneous  radar  cross  section  ak  of  realistic  targets  strongly  depends  on  the  radar  frequency  used 
and  the  current  aspect  angle.  For  this  reason,  statistical  models  are  used  for  describing  the  backscattering 
properties  of  the  targets.  In  many  practical  cases  is  described  by  gamma-densities: 

p((7k\d,m)  =  Gm(ak\  d) 

=  W*)"'  _»-t  _-«*(»/*) 

T(m)  k 

In  this  equation  d  denotes  the  mean  RCS  of  the  target  that  is  usually  unknown,  but  constant  in  time  and  char¬ 
acteristic  of  a  certain  class  of  targets.  The  individual  samples  <t/(  are  assumed  to  be  statistically  independent 


(1) 

(2) 
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for  subsequent  dwells  (guaranteed  by  frequency  decoiTelation,  e.g.).  The  cases  m  =  1,2  are  referred  to  as 
Swerling-I  and  -III  fluctuations [3]. 

Let  the  instantaneous  target  signal  v*  =  (vi,  V2)  be  additively  corrupted  by  GAUSSian  noise  with  variance 
(7^.  As  the  signal  components  are  assumed  to  be  statistically  independent,  the  pdf  of  the  resultant  signal 
s*  =  (.s  i ,  S2)  is  therefore  given  by  a  product  of  GAUSSians: 

p(SkWk)  =  1;  vi,  al)  JsT{sr,  v2,  a2).  (3) 

The  normalized  scalar  quantity  a2  =  (s2+s2) /2o2  derived  from  the  vector  signal  s*  is  RlCEian  distributed [3]: 

p(a2k  |snfe)  =  e~ak~snk  I0[2akVsr^)  with:  sn*  =  (vj  +  v11)l2c1n.  (4) 

Hence,  sn*  denotes  the  instantaneous  signal-to-noise  ratio  of  the  target  being  proportional  to  ak .  Its  expec¬ 
tation  with  of  respect  to  p(a2k\snk)  is  E[ajr]  =  1  +  sn*;  pure  noise  (sn*  =  0)  has  thus  unit  power.  Due  to  the 
RCS  model  previously  discussed,  sn*  is  gamma-distributed  with  the  mean  SN:  p(sn*  |SN)  =  <5m(sn*;  SN). 
The  conditional  density  of  a2k  given  SN  is  thus  obtained  by  [1]: 

p(a2k |SN)  =  {“dsn*  p(a2k\snk)  p(sn*|SN)  (5) 

_  /  m+SN  \  ~ m  „-mal/(m+ SNt  r  (  ~alSN  \ 

~\  m  )  e  (  m+SN  )• 

where  Lm_  1  denotes  the  Lagueme  polynomials.  For  Swerling-I/III:  Lq(—x)  =  1,  L\(—x)  =  1  +  x.  Evidently, 
p(ak  |SN)  can  be  interpreted  as  a  gamma  mixture  with  E[ a k  \  =  1  +  SN. 

1.2.3  Signal-to-Noise  Ratio 

Any  sensor  model  for  phased-array  tracking  has  to  provide  a  functional  relationship  between  the  expected 
signal-to-noise  ratio  SN*  at  time  tk,  the  sensor  parameters  (transmitted  energy,  beam  position)  and  the  target 
parameters  (mean  RCS,  target  position).  With  a  GAUSSian  beam  form[8]  and  using  the  radar  range  equation 
we  assume: 


SN,  =  SN„  (i)  (ft)  (a)-4  e-2"*  (7) 

with  Abit  =  (Uk  -  buk)2 / Bk  +  (Vft  -  bvk)2/B2  (8) 

r*  is  the  actual  target  range  at  time  r*  while  uk,  v*  denote  the  related  direction  cosines.  With  the  beam 
position  b*  =  (buk,  bvk)  and  the  (one-sided)  beam  width  B,  Ab*  is  a  measure  of  relative  beam  positioning 
error.  The  radar  parameter  SNo  is  the  expected  mean  signal-to-noise  ratio  of  a  target  with  a  standard  mean 
cross  section  do  at  a  reference  range  ro  that  is  directly  (Ab*  =  o)  illuminated  by  the  beam  with  the  energy 
Co-  Due  to  the  functional  relationship  stated  in  Equation  6  the  signal  strength  a2k  can  be  interpreted  as  a 
measurement  of  a. 

1.2.4  Detection/Measurements 

A  detection  is  assumed  if  the  received  signal  strength  exceeds  a  certain  detection  threshold:  a2k  >  Ad-  For  a 
given  m  in  the  fluctuation  model  (Equation  2),  the  detection  probability  Pd  is  a  function  of  SN  and  Ad'- 

Pd( SN,  Ad,  m)  =  \^da\  p{a\ |SN).  (9) 

The  false  alarm  probability  Pp  is  analogously  obtained: 


Pf(Ad)  =  Pd( 0,  AD,m)  =  eTiB. 


(10) 
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Integration  results  in  explicit  expressions  for  Pd  [3].  For  Swerling-I/III  fluctuations  we  obtain: 

P'(SN,  Ad)  =  e“TTM  =  (11) 

P"/(SN,AD)  =  e-*(l  +  |g^).  (12) 

9  9 

For  target  tracking  ak  is  available  after  a  detection,  i.e.  ak  >  Ad-  We  thus  need  the  conditional  density: 


p{a2k\a\  >  Ad.SN ,m)  = 


( p(al\SN)/PD(SN,  Ad,  m) 

1° 


for  a~k  >  Xd 
for  a2k  <  Ad 


(13) 


For  strong  targets 


(see  Equation  6). 
Pd(SN,  Ad,  m)  ~ 


with 


we  can  assume  SN  «  1  +  SN 


m  +  SN  and  thus  approximately  obtain: 


P(a2k  |SN) 


(|)-e^Vl(-«{) 


(14) 


On  the  other  hand,  let  the  detection  probability  for  m  ^  1  be  approximately  given  by: 
P^(SN,  Ad)  (i-e.  Swerling  I).  We  can  therefore  write: 


Sm(crk,  SN)  = 


p(a\a2  >  Ad,  SN,  tn)  a  Sm(a2;  SN) 

_  /  (f )~"  e-(^+^VSN  Lm_x  (  _  for  a2  >  Ad 

)  for  a?  <  Ac. 


(15) 

(16) 


Let  us  assume  that  monopulse  localization  after  detection  result  in  bias-free  measurements  u'k,  v'k  of  the 
direction  cosines  and  range  with  GAUSSian  measurement  errors.  According  to  [3]  the  standard  deviations 
(juk’v  depend  on  the  beam  width  B  and  the  instantaneous  sru  in  the  following  manner:  tr"  '  oc  B / ^/sn /<  a 

P/V ' a2k  —  1.  As  sn/<  is  unknown,  in  the  last  approximation  a2k  is  used  as  a  bias-free  estimate  of  sn/<  (E[rq  ]  = 
1  +  sn^).  The  range  error  is  assumed  to  be  GAUSSian  with  a  constant  standard  deviation  ur .  Evidently,  this 
model  of  the  measurement  process  does  not  depend  on  the  RCS  fluctuation  model. 


1.2.5  IMM  Dynamics  Model 

Systems  with  MARKOVian  switching  coefficients  [2]  are  well-suited  for  modeling  the  different  phases  in  a 
mission.  They  are  represented  by  multiple  dynamics  models  with  a  given  probability  of  switching  between 
the  models.  Let  4  denote  the  dynamics  model  indexed  by  4  at  time  4 .  The  transition  probabilities  are  thus 
part  of  the  modeling  assumptions,  however  a  priori  unknown  in  a  real  application.  Fortunately,  the  tracking 
performance  does  not  seem  to  depend  critically  on  their  particular  choice  [4].  The  related  Markov  process 
is  described  by: 

Pfak,  ikl^k— 1  ■  ffc— l)  —  Pikik- 1  W"(x^,  Q  jk),  1  <  i  k ,  4  — !  5-i  r.  (12) 

Hence,  the  model  jump  process  is  a  Markov  chain  with  model  transition  probabilities  pikik =  p(4|4-i), 
the  individual  models  being  linear-GAUSSian.  F(/t  denotes  the  evolution  matrix,  while  Di(.  reflects  the  ma¬ 
neuvering  capability  of  the  model  indexed  by  4  • 

1.3  Adaptive  Track  Maintenance 

As  discussed  in  the  previous  talk,  tracking  is  an  iterative  updating  scheme  for  conditional  probability  den¬ 
sities  p(Xk\Zk)  that  describe  the  current  target  state  x^  given  all  available  resource  allocations  Zk  and  the 
underlying  a  priori  information  in  terms  of  statistical  models.  Essentially,  the  update  consists  of  a  prediction 
step  followed  by  filtering.  The  knowledge  of  the  target  state  at  time  4  before  a  new  allocation  has  taken 
place  is  given  by  p(Xk+ 1  \Zk)-  Allocation  decisions  for  a  certain  time  must  thus  be  based  on  this  density  that 
essentially  depends  on  the  underlying  dynamics  model. 
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1.3.1  BAYESian  IMM-MHT 

Due  to  the  Total  Probability  Theorem,  dynamics  modeling  according  to  IMM  results  in  finite  mixture  densi¬ 
ties: 


p{xk\Zk)  =  2  Mik.  Zk )  p{h\Zk)  (18) 

h 

with  ik  =  {ik, ...  ,i  i}  denoting  a  particular  model  history,  i.e.  a  sequence  of  possible  hypotheses  regard¬ 
ing  the  target  dynamics  model  from  the  initial  observation  up  to  the  most  recent  measurement  at  time  t k. 
This  is  in  analogy  to  the  discussion  in  the  previous  talk.  As  the  number  of  terms  in  Equation  1 8  exponen¬ 
tially  increases  with  increasing  k,  various  techniques  have  been  developed  for  approximately  representing 
p(xk\Zk)  by  mixtures  with  a  constant  number  of  components  at  each  tk.  With  i*  =  {ik, . . . ,  ik-n+ 1 }  denoting 
a  sequence  of  possible  model  hypotheses  “ n  scans  back”,  we  are  looking  for  approximations  by  Normal 
Mixtures, 

P(xk\Zk)  ~  H  XjJ,  Pjj),  (19) 

•k 

1/7 

with  pik  =  p(ik\Zk)  and  J*f(xkm,  Xj*,  Pj*)  =  p(xk\ik,  Zk).  Due  to  Bayes’  Rule,  x^  and  P^  of  each  mixture 
component  are  iteratively  obtained  by  formulae  essentially  based  on  Kalman  filtering.  Also  the  weighting 
factors  p^  obey  simple  formulae[2,  16].  In  case  of  a  single  dynamics  model  (r  =  1),  the  conditional  densities 
p{xk\Zk)  are  strictly  given  by  GAUSSians;  i.e.  the  filtering  loop  is  simply  KALMAN  filtering. 

For  n  =  1,  p(xk\Zk)  is  approximated  by  a  mixture  with  r  components  according  to  the  r  dynamics 
models  used.  GPB2  and  standard  IMM  algorithms  are  possible  realizations  of  this  scheme  [2].  For  standard 
IMM  the  approximations  are  made  after  the  prediction,  but  before  the  filtering  step,  while  for  GPB2  they 
are  applied  after  the  filtering  step.  Hence,  GPB2  requires  some  more  computational  effort.  More  generally 
speaking,  the  basic  idea  of  standard  IMM  (moment  matching  directly  after  the  prediction  step)  may  easily 
be  adopted  to  n  >  2  providing  a  more  accurate  approximation  of  the  densities  p{xk\Zk)  than  for  n  =  2.  At 
each  stage  k  of  the  tracking  loop  the  parameters  pik,  x^,  approximately  represent  the  density  p(xk\Zk). 


1.3.2  RCS  as  a  State  Variable 

Fet  us  treat  the  normalized  mean  RCS  of  the  target,  sk  =  ak/&o,  as  an  additional  component  of  the  state 
vector.  As  the  signal  strength  after  detection  may  be  viewed  as  a  measurement  of  sk,  let  us  consider  the 
conditional  density  p{xk,  sk\Zk)  =  p(sk\xk,  Zk)p{xk\Zk)-  The  calculation  of  p{xk\Zk)  was  discussed  above. 
For  p(sk\xk,  Zk)  Bayes’  Rule  yields  up  to  a  normalizing  constant: 

p(sk\xk,a2k,Zk-l)*Sm(a2k-  SN)  p{sk\xk,  Zk~l).  (20) 

Fet  us  assume  that  p{sk\xk,  Zk~l)  are  given  by  inverse  gamma  densities: 

p(sk\xk,Zk~l)  =  (s*;  s*|*-i)  (21) 

with:  s)  =  [((//  -  l)Sy7T0/)]  (22) 

defined  by  s  =  E[  .v  |  >  0  and  a  parameter  p  >  1.  For  p  >  2  the  related  variance  exists:  V[s]  =  s2/(p  -  2). 
This  class  of  densities  is  invariant  under  the  successive  application  of  Bayes  Rule  according  to  Equation  20, 
for  up  to  normalization  we  obtain: 


^ in  (  ^4  ’  ^19)  |t_i  ( Sk ,  l)  OC 


— ttl  f^k\k—l  ktl  1 

ak  Sk 


exp( 


(fik\k-l~l)h\k-l  + 
Sk 


7 

mak+A  ft 


■)  a  xnk(sk\  h) 


(23) 
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with: 


=  S  No(t)(^rV2^ 

/v  _  Mk\k-\  —  1  /v  ,  (mak+Aj))/ ait 

Sk  ~  lik\k-\+m-\  Sfclfc_1  /(/t|fc-l+'«— 1 

Pk  =  Bk\k-\  +  m- 


(24) 

(25) 

(26) 


With  reference  to  s k  the  density  lflk  (sk ;  Sk)  is  correctly  normalized.  Due  to  ak  =  ak{rk,  uk ,  v/<),  however, 
it  depends  on  the  target  position.  In  order  to  preserve  the  factorization  of  p(\k ,  sk\Zk)  in  a  normal  mixture 
related  to  x*  and  an  inverse  gamma  density  related  to  sk  we  use  the  approximation: 

ak  »  (a)  (ii)-4  e-2f(^-6p2+(v,-6p2}/s2  (27) 

ryt,  Ufc,  v,t  are  the  MMSE  estimates  for  rk,  uk  and  v*  derived  from  p(x*  |Zfc).  Hence,  ak  compensates  both,  the 
estimated  positioning  error  of  the  radar  beam  and  the  propagation  loss  due  to  the  radar  equation.  Assuming 
sk  to  be  constant,  we  have  IMk]k  (,s7;  sk\k- 1 )  =  ^nk\k-i  isk’  $k- 1  )•  In  principle,  a  dynamics  regarding  the  cross 
section  might  be  introduced. 


1.4  Adaptive  Sensor  Control 

Based  on  the  previous  considerations  adaptive  techniques  for  combined  tracking  and  sensor  management  are 
discussed,  i.e.  control  of  data  innovation  intervals,  radar  beam  positioning,  and  transmitted  energy  [17,  19]. 


1.4.1  Revisit  Time  Control 

The  time  tk  of  a  radar  allocation  Rk  is  determined  by  the  current  lack  of  information  conveniently  described 
by  the  covariance  matrix  P*!*— i  of  the  predicted  state  estimate  \k\k-i  [8].  Since  the  predicted  pdf, 

p(Xk\Zk-l)  =  ^  Rk\k-1  •^'(Xfc’  XAr|/c-l’  1 ^8) 

;  k 


is  a  normal  mixture,  we  obtain: 

^  Rk\k-1  (P/t|fc-l  +  (XA:|/t-l 


xk\k-l)  (xjt|Jfc— 1 


l  )T )  • 


(29) 


The  covariance  of  the  individual  mixture  components  the  faster  grow  in  time  the  more  often  maneuvers  were 
assumed  in  the  corresponding  model  histories.  This  has  impact  on  the  total  covariance  P*:|fc— 1  according  to 
the  corresponding  weighting  factors.  In  addition  P*|*— 1  is  “broadened”  by  the  positive  definite  spread  term. 
Evidently,  the  adaptive  IMM  modeling  affects  P*|fc_i  in  a  rather  complicated  way.  A  scalar  measure  of  the 
information  deficit  is  provided  by  the  largest  eigenvalue  of  the  covariance  of  the  predicted  target  direction 
(in  terms  of  u,  v).  Let  it  be  denoted  by  Gk\k-\-  A  track  update  is  allocated  when  the  Gk\k- 1  exceeds  a 
predetermined  proportion  of  the  squared  radar  beam  width  B  : 

Gk\k-i  >  (v0B)2.  (30) 

The  relative  track  accuracy  vo  introduced  by  this  criterion  is  a  measure  of  the  minimum  track  quality  required 
and  a  parameter  to  be  optimized.  In  many  practical  applications  vq  =  0.3  is  a  reasonable  choice  [8]. 


1.4.2  Transmitted  Energy 

In  view  of  the  tracking  system  the  sensor  performance  is  mainly  characterized  by  the  signal-to-noise  ratio  that 
determines  both,  the  detection  probability  and  the  measurement  error.  By  suitably  choosing  the  transmitted 
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energy  per  dwell  ek,  the  expected  signal-to-noise  ratio  SN^_i  can  be  kept  constant  during  tracking.  Besides 
vo,  SN*;|k_i  is  an  additional  parameter  subject  to  optimization.  Since  vo  may  be  viewed  as  a  measure  of  the 
beam  positioning  error,  the  energy  ek  at  time  tk  is  defined  by  this  condition  (Equation  8): 

SNfc|fc_i  =  const.  (31) 

-  (32) 

By  this  particular  choice  the  influence  of  the  radar  range  equation  is  compensated  (at  least  for  a  certain 
range  interval).  For  the  mean  radar  cross  section  a  either  a  worst-case  assumption  or  estimates  from  target 
amplitude  information  can  be  used.  Also  the  track  quality  vo  affects  the  transmitted  energy.  As  a  side  effect 
of  this  choice,  the  standard  deviations  oy  ’v  of  the  u,  v -measurements  arc  kept  constant  on  an  average. 

1.4.3  BAYESian  Local  Search. 

Intelligent  algorithms  for  beam  positioning  and  local  search  are  crucial  for  IMM-type  tracking.  Too  simple 
strategies  may  easily  destroy  the  benefits  of  the  adaptive  dynamics  model,  because  track  loss  immediately 
after  a  model  switch  is  to  be  expected.  For  this  reason  the  optimal  approach  based  on  the  predicted  densities 
p(x-k\Zk~l)  proposed  in  [8]  is  adopted  to  IMM  tracking.  The  beam  position  of  the  first  dwell  is  simply 
given  by  predicted  direction  u^-i,  v*|*_i  derived  from  p(x/<  ).  If  no  detection  occurs  in  the  first  dwell, 

however,  this  very  result  provides  useful  information  on  the  target  state.  We  thus  have  to  calculate  the 
conditional  density  given  the  event  ->D  “no  detection  in  this  direction”.  Bayes’  Rule  yields: 

p{uk,vk\-V,Zk-l)*{\-PD{uk,vk))p(uk,vk\Zk-X)  (33) 

up  to  a  normalizing  factor.  In  this  density  the  detection  probability  Pd  depends  on  the  expected  signal- 
to-noise  ratio  (Equation  8).  By  this  the  fluctuation  model  for  the  target  RCS  enters  the  discussion  as  it 
determines  the  functional  forms  of  Pd  (Equation  12).  The  two  dimensional  density  p{uk,  vy  |->D,  Zk  *)  can 
easily  be  calculated  on  a  grid.  The  beam  position  for  the  next  dwell  is  then  simply  provided  by  the  maximum 
of  p{uk,  vk | - 1 D ,  Zk~l).  This  scheme  is  repeated  until  a  detection  occurs.  Since  the  maximum  is  searched,  the 
computation  of  the  normalization  integral  is  not  required.  Numerically  efficient  realizations  are  possible. 

Figure  2  illustrates  this  scheme  for  one  of  the  examples  discussed  below.  First  the  predicted  target  pdf  (a 
mixture)  is  shown  as  a  function  of  its  direction  coordinates.  With  high  probability  the  target  is  expected  to  be 
in  the  dark  region,  the  true  target  position  being  indicated  by  a  dark  grey  dot.  The  light  dot  denotes  the  beam 
position  of  the  next  dwell.  The  related  detection  probability  is  26%.  However,  no  detection  occurred  during 
the  first  dwell.  We  thus  calculate  the  conditional  pdf  given  that  event.  The  resulting  pdf  differs  significantly. 
The  previous  maximum  decreased  in  height,  while  the  global  maximum  is  at  a  different  location.  Again 
no  detection  occurred;  the  resulting  pdf  is  showed  next.  Now  the  algorithm  decides  to  look  again  near  the 
position  at  dwell  1 .  In  addition,  two  smaller  local  maxima  appear  that  increase  in  size  at  the  next  dwell  as  no 
detection  occurred.  We  finally  obtain  a  decision  which  leads  to  success.  The  last  picture  shows  the  updated 
pdf. 

1.5  Discussion:  Simulated  Examples 

Simulation  results  provide  hints  to  what  extend  the  total  performance  of  multiple-target  air  surveillance  by 
phased-array  might  be  improved  by  using  adaptive  techniques  for  combined  tracking  and  sensor  control. 
Four  questions  are  addressed: 

1.  Which  resource  savings  (allocation  time,  energy)  can  be  expected  by  using  adaptive  dynamics  models? 

2.  How  should  the  IMM  dynamics  modeling  be  designed  (e.g.  number  of  models,  transition  matrix)? 
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Figure  2:  BAYESian  local  search:  five  consecutive  dwells 
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Figure  3:  Horizontal  projections  and  kinematical  quantities 


3.  Which  energy  savings  can  be  expected  by  exploiting  target  amplitude  information  for  sensor  control? 

4.  Why  is  BAYESian  local  search  important  when  IMM  dynamics  models  are  used  for  revisit  time  con¬ 
trol? 


In  general  we  follow  the  parameter  and  threshold  settings  recommended  in[8] .  To  exclude  false  alarms 
due  to  receiver  noise,  the  false  alarm  probability  is  Pp  =  I  (F1 .  False  returns  due  to  clutter  or  ECM  are  not 
considered.  The  error  of  range  measurements  (standard  deviation)  is  oy  =  100  m,  the  radar  beam  width 
B  =  1°.  We  assume  a  minimum  time  interval  of  20  ms  between  consecutive  dwells  on  a  particular-  target 
and  statistically  independent  signal  amplitudes  (frequency  decorrelation,  e.g.).  The  reference  range  is  set 
to  ro  =  80  km.  Antenna  coordinates  (direction  cosines,  range)  are  used  also  for  tracking;  non-linearities 
introduced  by  these  no-Cartesian  coordinates  are  taken  into  account  [7].  The  maneuvering  capability  of  the 
targets  is  thus  characterized  by  two  parameters:  maneuver  correlation  time  8  and  acceleration  width  Z.  For 
r  =  2,  3  we  consider  the  parameter  sets: 
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Mi 

(worst-case  model): 

Ei  =  60  m/s2, 

9\  =  30  s 

Mi 

(best-case  model): 

Z2  =  1  m/s2, 

92  =  10  s 

(Pij)  ~  (  2  .9  )  ’  (Pij') 

m3 

(medium-case  model): 

X3  =  30  m/s2, 

02  =  30  s 

We  observed  that  the  performance  does  not  critically  depend  on  the  particularly  switching  probabilities  p(/ 
chosen.  A  detailed  mismatch  analysis,  however,  has  not  been  performed.  A  track  is  considered  to  be  lost  if 
more  than  50  dwells  occur  in  the  local  search  of  if  the  beam  positioning  error  A bk  is  greater  than  3 B.  We 
thus  permit  even  a  rather  extensive  local  search  that  correspondingly  burdens  the  total  energy  budget.  In  all 
simulations  considered  below  (1000  runs)  the  relative  frequency  of  track  loss  is  less  than  2%. 

1.5.1  Benchmark  Trajectories 

The  horizontal  projection  of  four  standard  benchmark  trajectories  (military  cargo  aircraft,  medium  bomber, 
fighter/attack  aircraft,  and  anti-ship  missile)  is  shown  in  Figure  3  along  with  representative  kinematical  char¬ 
acteristics  such  as  acceleration  (solid  line),  range  (dashed),  height  (dotted),  and  speed  (solid).  They  have 
been  proposed  in  [5,  6]  and  cover  a  rather  wide  range  of  militarily  relevant  targets.  The  missile  trajectory 
might  serve  to  explore  the  performance  limits  of  the  algorithms.  In  principle,  missiles  can  execute  even 
stronger  maneuvers.  It  is  questionable,  however,  if  for  those  objects  and  their  individual  missions  the  dy¬ 
namics  models  discussed  above  remain  applicable.  All  targets  are  tracked  over  a  period  of  180  s.  The  RCS 
fluctuations  are  described  by  a  Swerling-III  model.  The  mean  cross  sections  significantly  vary  from  target  to 
target  (4.,  2.,  1.2,  .5  nr). 

The  discussion  is  confined  to  a  few  intuitively  clear  and  simple  performance  measures  obtained  by 
Monte-Carlo  simulation  (1000  runs).  In  general  a  single  performance  measure  is  not  sufficient  as  there 
may  exist  applications  where  the  transmitted  energy  is  the  limiting  factor,  while  in  a  different  scenario  the 
number  of  radar  allocations  must  be  kept  low. 

The  adaptivity  becomes  visible  if  the  performance  is  evaluated  as  a  function  of  the  tracking  time  that  can 
be  compared  with  the  kinematics  of  the  individual  trajectories  (Figure  3).  Here  we  used  histograms  with  100 
cells.  In  particular-  we  considered:  the  mean  revisit  intervals,  the  mean  number  of  dwells  for  a  successful 
update,  the  mean  number  of  sensor  allocations  totally  required  for  track  maintenance,  the  mean  energy  spent 
for  a  successful  allocation,  the  mean  energy  totally  spent  for  track  maintenance,  and  the  mean  RCS  of  the 
targets  estimated  during  tracking. 

Four  tracking  filters  were  compared:  worst-case  Kalman  filter  (KF),  standard  IMM  filter  with  two  or 
three  models,  respectively  (S-IMM2  3),  and  IMM-MHT  filtering  with  model  histories  of  length  n  =  4.  For 
IMM-MHT  with  n  >  4  the  performance  characteristics  only  slightly  change.  We  thus  conclude  that  n  =  4 
already  provides  a  good  approximation  to  optimal  filtering  (at  least  for  the  scenarios  considered  here).  With 
reference  to  target  amplitude  information  we  considered  three  cases:  1)  the  target  RCS  a  is  known  and 
used  for  energy  management.  2)  The  mean  RCS  a  is  unknown  and  to  be  estimated  during  tracking.  3)  A 
worst-case  assumption  is  used  for  all  targets  (a  -  0.5  m2). 

1.5.2  Design  of  IMM  modeling 

Practically,  the  question  arises  how  many  models  should  be  used  in  the  IMM  approach.  In  addition  it  must  be 
clarified  if  each  trajectory  needs  an  individual  modeling  or  if  the  same  IMM  modeling  can  be  used  without 
significant  loss  of  performance.  For  the  fighter  scenario  a  worst/best-case  model  should  be  appropriate  at  first 
sight.  Trajectory  1  (Cargo  Aircraft),  however,  shows  that  military  targets  can  occur  for  which  medium-case 
models  are  sufficient.  To  answer  these  questions  we  used  IMM  with  two  (r  =  2,  Mi,  Mi)  and  three  models 
(r  =  3,  Mi,  Mi,  M3),  respectively,  with  vo  =  0.3,  SN^—i  =  50.  In  which  way  these  parameters  affect  the 
performance  is  discussed  further  below.  Figure  4  (first  row)  shows  the  resulting  mean  revisit  intervals  for  all 
trajectories.  The  kinematical  target  characteristics  are  clearly  mirrored.  We  observed: 
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Figure  4:  Revisit  intervals,  allocations,  and  energy  for  different  filters  and  trajectories 
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1.  As  expected,  Kalman  filtering  (r  =  1,  Mi)  leads  to  constant  revisit  intervals  that  are  comparable  for 
all  trajectories.  This  is  no  longer  true  for  S-IMM.  The  resultant  curves  related  to  r  =  2  (solid)  and 
r  =  3  (dashed)  significantly  differ  from  each  other.  The  onset  of  maneuvers  (Figure  3)  strongly  affects 
the  mean  update  intervals  and  thus  illustrates  the  adaptivity  of  the  algorithm. 

2.  The  difference  between  the  cases  r  =  2,  3  vanishes  however,  if  IMM-MHT  is  used.  If  model  histories 
are  permitted  (here  n  =  4),  it  seems  to  be  irrelevant  if  besides  worst/best-case  assumption  additional 
medium-case  models  are  used.  Even  longer  histories  or  further  models  (r  >  3)  do  not  significantly 
improve  the  performance  obtained  with  r  =  2  and  n  —  4.  For  a  suitable  (!)  choice  of  the  switching 
probabilities  the  performance  of  S-IMM4  approaches  closes  to  B-EVOVU;  for  B-IMM4  no  improvement 
over  B-IMM2  was  observed. 

3.  For  the  bomber  and  the  fighter  S-IMM3  (Mi ,  M2,  M3)  outperforms  S-IMM2  (Mi ,  M2),  in  spite  of  the 
fact  that  for  these  trajectories  worst-case  maneuvers  occur  only  and  the  medium-case  model  appealed 
to  be  unnecessary  at  first  sight.  The  difference  between  r  =  2  and  r  =  3,  however,  is  not  as  clear  as  for 
scenario  1  (cargo  aircraft). 

4.  For  the  moderately  maneuvering  cargo  aircraft  the  question  arises  if  the  performance  can  be  im¬ 
proved  by  using  a  medium/best-case  IMM  modeling.  We  found  that  worst/best-case  IMM-MHT  and 
medium/best-case  IMM-MHT  differ,  but  not  very  much.  This  indicates  that  worst/best-case  IMM- 
MHT  has  a  more  or  less  “universal”  character,  i.e.  it  does  not  critically  depend  of  the  scenario 
considered  (at  least  within  certain  limits). 

These  observations  indicate  that  the  mixtures  p(xk\Zk)  for  n  =  4,  r  =  2  have  enough  internal  degrees  of 
freedom  to  provide  an  adequate  representation  of  the  actual  target  behavior.  Refined  approximations  by 
even  more  mixture  components  seem  to  be  irrelevant  for  the  trajectories  considered.  A  Rule  of  thumb:  A 
worst/best-case  analysis  of  the  problem  along  with  IMM-MHT  seems  to  be  sufficient  to  achieve  a  nearly 
optimal  tracking  performance.  Evidently,  for  two  dynamics  models  reasonable  and  intuitive  assumptions 
for  the  switching  probabilities  are  easily  obtained.  IMM-MHT  thus  enables  a  more  simplified  dynamics 
modeling  without  significant  loss  of  performance. 

1.5.3  Gain  by  IMM  modeling 

For  investigating  the  gain  by  adaptive  dynamics  models,  let  us  for  the  present  assume  that  the  mean  RCS 
of  the  target  is  known  and  used  for  energy  management.  Figure  4  shows  the  mean  number  of  allocations 
required  for  track  maintenance  (KF,  S-IMM2,  B-IMM2).  As  expected,  for  KF  the  mean  number  of  revisits 
linearly  increases  with  increasing  tracking  time  and  is  nearly  the  same  for  all  trajectories.  By  adaptive 
dynamics  modeling,  however,  the  number  of  sensor  allocations  is  reduced. 

1.  Compared  with  KF,  IMM  results  in  significant  resource  savings.  There  is  an  improvement  by  IMM- 
MHT  over  S-IMM;  the  difference,  however,  is  less  significant  than  between  S-IMM  and  KF.  Besides 
simplified  modeling  assumptions,  the  practical  use  of  IMM-MHT  therefore  consists  in  the  exploration 
of  the  limiting  bounds  for  performance  improvements. 

2.  The  largest  gain  is  observed  for  the  cargo  aircraft  and  the  bomber.  In  case  of  the  fighter  the  allocations 
required  arc  reduced  by  about  50  %  compared  with  worst-case  Kalman  filtering.  Even  during  the  7 
g  weaving  of  the  missile  some  advantages  of  the  IMM  modeling  can  be  observed. 

Figure  4  (second  row)  also  shows  the  mean  number  of  dwells  per  revisit.  Up  to  peaks  corresponding  with 
the  onset  of  maneuvers,  it  is  constant  and  roughly  equal  for  all  filters  and  trajectories.  The  peaks  are  the 
higher  the  more  adaptivity  the  filter  shows,  i.e.  the  larger  the  revisit  intervals  can  be  during  inertial  flight. 
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Figure  5:  On  the  quality  of  RCS  estimates 


The  peaks  thus  indicate  that  for  abrupt  maneuvers  a  local  search  might  be  required.  This  is  the  price  to  paid 
for  increased  adaptivity.  Evidently,  intelligent  algorithms  for  beam  positioning  and  local  search  are  essential 
for  IMM  phased-array  tracking. 

These  observations  are  consistent  with  Figure  4  (third  row)  which  shows  the  mean  energy  spent  for 
track  maintenance  (relative  units).  Besides  the  target  maneuvers  these  curves  are  influenced  also  by  the 
current  target  range  (Figure  3,  dotted  line).  In  addition  the  mean  energy  spent  per  revisit  is  displayed.  Up  to 
characteristic  peaks  the  energy  per  revisit  is  roughly  the  same  for  all  tracking  filters. 

1.5.4  Quality  of  RCS  Estimates 

In  a  practical  application  the  mean  RCS  of  the  targets  to  be  tracked  is  unknown  and  might  be  estimated  from 
target  amplitude  information.  In  general,  the  estimators  used  should  be  at  least  approximately  bias-free,  the 
estimated  error  and  the  empirical  error  should  be  roughly  identical,  and  the  estimators  should  show  a  certain 
robustness  against  model  mismatch.  As  indicated  by  Figure  5,  the  estimator  previously  proposed  provides 
rather  satisfying  results  for  all  trajectories.  Using  IMM-MHT  for  tracking,  the  recursion  was  initiated  with 
a  =  .5  m2  (worst-case  assumption)  and  mo  =  1.01. 

The  solid  lines  show  the  mean  RCS  estimates  as  a  function  of  the  tracking  time.  For  all  scenarios  it  is 
roughly  constant  and  corresponds  with  the  actual  values  (4,  2,  1.2,  0.5  m2).  The  dotted  lines  indicate  the 
mean  estimation  error  (available  in  the  simulation).  The  curves  show  peaks  that  are  related  to  the  onset 
of  maneuvers  and  the  corresponding  lack  of  track  accuracy.  The  dashed  lines  denote  the  mean  standard 
deviation  calculated  by  the  estimator  itself.  Tracking  and  RCS  estimation  are  closely  interrelated:  Only 
when  tracked  over  a  certain  period  of  time,  estimates  are  reliable  enough  to  distinguish  between  the  target 
classes.  A  satisfying  RCS  estimation  by  signal  processing  only,  i.e.  without  a  temporal  integration  along  the 
estimated  trajectory  does  not  seem  to  be  possible.  In  this  context  IMM  retrodiction  techniques  [16]  might 
be  considered  which  can  provide  more  accurate  estimates  of  the  trajectory  and  thus  more  accurate  RCS 
estimates. 

1.5.5  RCS  Model  Mismatch 

The  backscattering  properties  of  real  targets  arc  highly  complex.  A  practicable  method  for  estimating  the 
RCS  must  thus  show  some  robustness  against  model  mismatch.  To  get  a  first  hint,  we  generated  in  our 
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target  type 

processed 

simulated 

RCS  [m2] 

estimated  error 

estimation  error 

energy 

Bomber 

III 

III 

1.96 

0.22 

0.28 

0.37 

III 

I 

2.34 

0.25 

0.57 

0.42 

I 

III 

1.77 

0.28 

0.33 

0.39 

I 

I 

2.03 

0.32 

0.34 

0.44 

Fighter 

III 

III 

1.19 

0.13 

0.16 

0.82 

III 

I 

1.41 

0.15 

0.34 

1. 

I 

III 

1.07 

0.16 

0.19 

0.86 

I 

I 

1.22 

0.18 

0.19 

0.99 

Table  1 :  Mismatch  regarding  the  fluctuation  model 


simulation  amplitude  information  according  to  both  Swerling  I  and  III  being  processed  according  to  both 
modeling  assumptions.  The  results  for  the  four  possible  combinations  are  summarized  in  table  1.  Besides 
the  quantities  already  shown  in  Figure  5,  we  also  listed  the  total  energy  spent  for  tracking  (relative  units). 

1 .  For  matching  models  the  RCS  estimates  are  nearly  bias-free  and  more  or  less  roughly  consistent. 

2.  For  Swerling  III  fluctuations  the  estimates  are  more  accurate  than  in  case  of  Swerling  I. 

3.  For  Swerling  I  (no  mismatch)  more  energy  is  spent  than  for  Swerling  III  (keeping  SN^^_i  constant). 

4.  If  Swerling  I  amplitudes  are  processed  according  to  Swerling  III,  the  RCS  is  overestimated,  consis¬ 
tency  is  lost. 

5.  It  is  underestimated  if  Swerling  III  amplitudes  are  processed  according  to  Swerling  I. 

6.  Mismatch  does  not  greatly  affect  the  performance  (energy). 


1.5.6  Energy  Management 

Finally  we  have  to  show  up  to  what  degree  the  radar  energy  to  be  spent  can  be  reduced  by  estimating  the  RCS 
in  comparison  to  worst-case  assumptions.  In  Figure  5  the  mean  radar  energy  spent  for  track  maintenance 
is  displayed.  The  dotted  lines  refer  to  IMM-MHT  tracking  using  the  true  RCS  of  the  targets  (as  previously 
discussed).  In  a  practical  application  this  can  not  be  realized;  the  resultant  curves,  however,  may  serve  as 
a  reference  to  discuss  the  performance  of  RCS-adaptive  algorithms.  The  solid  lines  denote  methods  that 
exploit  signal  strength  information  for  estimating  the  RCS  (Worst-Case  Kalman  filter,  IMM-MHT).  By 
dashed  lines  algorithms  are  indicated  that  use  a  worst-case  assumption  (here:  0.5  nr,  missile)  on  the  RCS 
(KF,  IMM-MHT). 

A  comparison  between  sensor  control  by  using  the  true  RCS  (not  available  in  a  real  application)  and 
methods  exploiting  recursive  RCS  estimates  is  of  particular  interest.  The  largest  deviation  is  observed  for 
scenario  1  (d  =  4  m2).  This  is  to  be  expected,  as  the  recursion  was  started  with  a  worst-case  assumption. 
The  discrepancy  between  both  curves,  however,  is  not  very  significant  in  all  four  cases.  Compared  with 
IMM-MHT  (Worst-Case  RCS)  it  can  be  neglected.  The  difference  between  sensor  control  with  known 
and  estimated  RCS  is  roughly  constant  during  tracking.  We  thus  conclude  that  it  is  caused  primarily  in 
the  initiation  phase  where  not  much  signal  strength  information  is  yet  available.  As  soon  as  reliable  RCS 
estimates  have  been  produced,  the  performance  is  practically  identical.  Figure  5  also  shows  how  the  resource 
savings  due  to  adaptive  dynamics  models  and  RCS-adaptive  energy  management  are  related  to  each  other. 
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Figure  6:  Mean  Number  of  Allocations  for  Different  Filters  and  Trajectories 


target  type 

RCS 

filter 

AT  [s] 

revisits 

energie 

rel. 

Cargo 

worst  case 

Kalman 

3.3 

55.1 

2828 

8.5 

IMM-MHT 

6.5 

32.6 

1664 

5.0 

estimated 

Kalman 

2.0 

90.3 

750 

2.2 

IMM-MHT 

4.9 

45.0 

453 

1.4 

known 

4.9 

46.4 

334 

1 

Bomber 

worst  case 

Kalman 

3.1 

60.1 

5488 

5.5 

IMM-MHT 

6.4 

34.0 

2868 

2.9 

estimated 

Kalman 

2.2 

85.5 

2257 

2.3 

IMM-MHT 

5.1 

44.2 

1095 

1.1 

known 

5.0 

45.7 

993 

1 

Fighter 

worst  case 

Kalman 

2.8 

67.3 

5786 

2.6 

IMM-MHT 

4.9 

43.3 

3882 

1.7 

estimated 

Kalman 

2.2 

85.6 

3563 

1.6 

IMM-MHT 

3.9 

54.0 

2420 

1.1 

known 

3.8 

55.1 

2226 

1 

Missile 

worst  case 

Kalman 

2.1 

86.6 

8657 

1.3 

IMM-MHT 

4.1 

56.1 

6593 

1 

estimated 

Kalman 

2.2 

85.4 

9036 

1.4 

IMM-MHT 

4.2 

55.4 

7042 

1.1 

known 

4.1 

56.1 

6593 

1 

Table  2:  Gain  by  RCS-adaptive  energy  control 


In  table  2  for  all  scenarios  and  processing  methods  scalar  performance  measures  are  summarized:  target 
revisit  intervals  (AT),  sensor  allocations  required,  energy  spent  for  track  maintenance  (time  averages  taken 
over  the  tracking  time).  The  last  column  shows  the  energy  spent  by  the  various  methods  relative  to  IMM- 
MHT  with  known  RCS.  Compared  with  IMM-MHT  (Worst-Case  RCS)  the  gain  is:  3.8  (cargo  aircraft), 
2.6  (Bomber),  1.5  (Fighter),  .9  (anti-ship  missile).  Hence,  in  the  missile-scenario,  where  the  worst-case 
assumption  is  correct,  a  small  loss  of  performance  must  be  taken  into  account. 
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2  Sensor  Fusion  for  Ground  Target  Tracking 

Ground  surveillance  aims  at  near  real-time  production  of  a  dynamic  ground  picture.  This  task  comprises  track 
extraction  and  track  maintenance  of  single  ground  moving  vehicles  and  convoys,  mobile  weapon  systems 
or  military  equipment,  as  well  as  low-flying  targets  such  as  helicopters.  As  ground  target  tracking  is  a 
challenging  problem  all  available  information  sources  must  be  exploited,  i.e.  the  sensor  data  themselves  as 
well  as  background  knowledge  about  the  sensor  performance  and  the  underlying  scenario. 

For  long-range,  wide-area,  all-weather,  and  all-day  ground  surveillance  operating  at  high  data  update 
rates,  GMTI  radar  proves  to  be  the  sensor  system  of  choice  (GMTI:  Ground  Moving  Target  Indication). 
By  using  airborne  sensor  platforms  in  stand-off  ground  surveillance  applications  the  effect  of  topographical 
screening  is  alleviated,  thus  extending  the  sensors’  field  of  view.  In  [12]  characteristic  problems  of  signal 
processing  related  to  GMTI  tracking  with  STAP  radar  arc  discussed.  We  here  discuss  selected  tracking 
aspects,  which  arise  from  using  more  appropriate  sensor  models  and  expoiting  background  information.  The 
following  topics  are  of  particular  interest: 

•  Doppler-Blindness.  Even  after  platform  motion  compensation  by  using  STAP  techniques  [13],  ground 
moving  vehicles  can  well  be  masked  by  the  clutter  notch  of  the  sensor.  This  physical  phenomenon 
directly  results  from  the  low-Doppler  characteristics  of  ground  moving  vehicles  and  causes  interfering 
fading  effects  that  seriously  affect  track  accuracy  and  track  continuity.  Unless  appropriately  handled, 
Doppler-blindness  can  cause  serious  problems,  which  seem  to  be  even  more  difficult  in  the  presence 
of  Doppler  ambiguities. 

•  Road-Map  Information.  Even  military  targets  usually  move  on  road  networks,  whose  topographical 
coordinates  arc  known  in  many  cases.  Digitized  topographical  road  maps  such  as  provided  by  Geo¬ 
graphical  Information  Systems  (GIS)  should  therefore  enter  into  the  target  tracking  and  sensor  data 
fusion  process.  The  problem  of  battlefield  surveillance  is  not  considered  here  (i.e.  off-road  targets). 

•  Sensor  Data  Fusion.  Since  a  single  GMTI  sensor  on  a  moving  airborne  platform  can  record  the 
situation  of  interest  merely  over  short  periods  of  time,  sensor  data  fusion  proves  to  be  of  particular 
importance.  The  data  processing  and  fusion  algorithms  used  for  ground  surveillance  are  closely  related 
to  the  statistical,  logical,  and  combinatorial  methods  applied  to  air  surveillance. 

A  GMTI  radar  sensor  produces  measurements  of  kinematical  target  parameters  and  possibly  false  returns 
caused  by  residual  clutter,  e.g..  In  addition,  measurement  errors  and  sensor  parameters  are  provided.  The 
resulting  tracks  represent  the  kinematical  states  of  the  targets  along  the  corresponding  accuracies  and  as 
such  are  prerequisites  for  producing  a  ground  picture.  We  here  discuss  track  maintenance  for  well-separated 
vehicles  in  case  of  low  residual  clutter. 

The  papers  [9,  10]  arc  the  standard  references  explicitly  devoted  to  the  problem  of  GMTI  tracking  with 
single  and  multiple  sensors.  The  impact  of  sensor  modeling  on  GMTI  tracking  is  is  discussed  in  [18]. 

2.1  Discussion  of  an  idealized  Scenario 

Assuming  a  flat  earth  Figure  7  shows  an  idealized  scenario  with  two  airborne  GMTI  sensors  observing  a 
ground  vehicle  moving  with  constant  speed  (15  m/s  =  54  km/h)  parallel  to  the  x-axis  for  most  of  the  time. 
This  situation  is  typical  of  stand-off  or  gap-filling  ground  surveillance  missions.  In  the  second  half  of  the 
observation  period  over  Afmax  =  25  min  the  target  stops  for  7  min.  Then  it  speeds  up  again  reaching  its  initial 
velocity.  Finally,  the  target  leaves  the  field  of  view  of  sensor  2. 

In  Table  1  selected  sensor  and  platform  parameters  are  summarized.  hp,  vp  denote  the  constant  height 
and  speed  of  the  sensor  platforms  over  ground.  A r,  A cp  are  the  range  and  azimuth  regions  covered  by 
each  sensor  during  observation.  The  revisit  intervals  are  given  by  AT,  while  MDV  denotes  the  Minimum 
Detectable  Velocity,  a  GMTI-specific  sensor  parameter  being  important  to  ground  moving  target  tracking. 
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Sensor 

hp  [km] 

vp  [m/sec] 

A  r  [km] 

A  cp  [deg] 

AT  [sec] 

MDV  [m/sec] 

1 

10 

200 

[232,  292] 

[-128,  -67] 

15 

2 

2 

1 

40 

[22,  54] 

[  77,  172] 

10 

2 

Table  3:  selected  sensor  and  platform  parameters 


2.1.1  Doppler-Blindness 

Even  after  platform  motion  compensation  by  using  STAP-techniques  [13],  GMTI  tracking  may  remain  a 
difficult  task.  Unless  appropriately  handled,  in  particular  two  phenomena  can  cause  problems: 

1.  Sensor- to-target  geometries  can  occur  where  targets  to  be  tracked  are  masked  by  the  clutter  notch  of 
the  sensor.  This  results  in  a  series  of  missing  detections  until  the  geometry  is  changing  again. 

2.  As  stopping  targets  are  indistinguishable  from  ground  clutter,  the  early  detection  of  a  stopping  event 
itself  as  well  as  tracking  of  ‘stop  &  go’  targets  can  be  important  to  military  applications. 

The  impact  of  these  effects  on  the  detection  probability  is  shown  in  Figure  8  for  the  scenario  previously 
introduced.  For  both  sensors  we  observe  deep  notches  (dashed  line:  platform  1,  dotted  line:  platform  2). 
In  the  center  of  these  notches  the  radial  velocities  of  the  target  and  the  surrounding  ground  patch  are  very 
close  to  each  other,  thus  making  target  discrimination  by  Doppler  processing  (STAP  [13])  impossible.  This 
is  particularly  true  if  the  target  stops. 

The  dashed  and  solid  lines  in  Figure  9  denote  the  radial  velocities  of  ground  patches  around  the  target 
and  target  returns,  respectively.  The  area  shaded  in  grey  reflects  the  width  of  the  clutter  notchs  of  the  sensors, 
which  is  determined  by  the  individual  Minimum  Detectable  Velocities  (MDVs).  For  each  sensor  both  curves 
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Figure  8:  detection  probability 


are  closely  adjacent  to  each  other  indicating  that  the  target  is  moving  at  a  much  lower  speed  than  the  sensor 
platforms.  We  notice  sliding  intersections  between  the  curves.  They  are  responsible  for  the  relatively  long 
duration  of  Doppler-blind  phases. 

2.1.2  GMTI  Model 

For  dealing  with  this  phenomenon  we  propose  a  refined  model  for  describing  the  sensor  characteristics 
and  discuss  its  benefits  for  improving  the  performance  of  ground  target  tracking  and  sensor  data  fusion. 
The  model  is  adapted  to  STAP  techniques  in  that  the  detection  probability  assumed  in  the  tracking  process 
is  described  as  a  function  of  the  GMTI-specific  clutter  notch.  While  the  current  location  of  the  notch  is 
determined  by  the  kinematical  state  of  the  target  and  the  current  sensor-to-target  geometry,  its  width  is  given 
by  a  characteristic  sensor  parameter  (MDV). 

By  this  more  detailed  information  on  the  sensor  performance  can  be  incorporated  into  the  tracking  pro¬ 
cess.  This  in  particular  permits  a  more  appropriate  treatment  of  missing  detections.  In  other  words,  infor¬ 
mation  on  the  potential  reasons  that  might  have  caused  the  missing  detections  enters  into  the  tracking  filter. 
We  observed  that  by  this  measure  the  number  of  lost  tracks  can  significantly  be  reduced  while  the  track 
continuity  is  improved,  finally  leading  to  a  more  reliable  ground  picture. 


2.1.3  Sensor  Fusion 


Assuming  an  idealized  processing  architecture  (centralized  data  fusion),  the  mean  cumulative  revisit  interval 
A Tc  results  from  the  individual  revisit  intervals  AT),  i  =  1, . . . ,  ns,  of  ns  sensors  according  to 


1 

ATC 


1 


1 

AT,' 


(34) 
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Figure  9:  range  rate  (ground,  target) 


For  the  previous  example  (ns  =  2,  AT)  =  15  s,  AT)  =  10  s)  we  obtain  AT)  =  6  s.  The  mean  cumulative 
detection  probability  P‘D  referring  to  A Tc  is  given  by 


tV  ■  ATc 

(35) 
i=  1 

with  P'D  denoting  the  detection  probabilities  of  the  individual  sensors,  which  depend  on  the  corresponding 
sensor-to-target  geometries.  Evidently,  the  corresponding  revisit  intervals  AT)  also  enter  into  this  formula, 
which  describes  the  mean  improvement  of  the  overall  detection  performance  to  be  expected  by  sensor  data 
fusion.  The  larger  the  individual  revisit  interval  AT)  of  sensor  i,  the  smaller  is  the  effect  of  sensor  i  on  the 
collective  performance,  even  if  the  corresponding  individual  detection  probability  P‘D  is  large.  AT)  and  PLD 
are  averaged  quantities,  by  which  the  expected  performance  improvement  can  be  predicted  in  an  overall 
sense. 

Figure  9  shows  the  mean  cumulative  detection  probability  PCD  for  the  above  example  (solid  line).  The 
impact  of  the  clutter  notches  are  more  or  less  compensated.  Due  to  the  fact  that  PLD  is  related  to  the  mean 
cumulative  revisit  interval  ATC  =  6  s,  being  shorter  than  those  of  the  individual  sensors  (ATi  =  10  s,  AT)  = 
15  s),  Pld  is  smaller  than  the  detection  probability  of  the  sensor  dominating  at  that  time. 


2.2  Tracking  Preliminaries 

The  choice  of  a  suitable  coordinate  system  for  describing  the  underlying  sensor/target  geometry,  the  sensor 
platform  trajectory,  and  the  available  a  priori  information  on  the  dynamical  behavior  of  ground  moving 
targets  are  prerequisites  to  target  tracking. 
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2.2.1  Coordinate  Systems 

For  the  sake  of  simplicity  we  consider  three  coordinate  systems  in  which  the  underlying  physical  phenomena 
become  transparent: 

1.  Cartesian  ground  coordinates,  where  the  description  of  the  target  and  platform  kinematics  is  of  a 
particularly  simple  form, 

2.  the  moving  Cartesian  antenna  coordinate  system,  whose  x-axis  is  oriented  along  the  array  antenna  of 
the  GMTI  radar  mounted  on  the  airborne  sensor  platform, 

3.  the  sensor  coordinate  system,  in  which  the  measurements  of  the  kinematical  target  parameters  are 
described  (target  range,  azimuth,  and  range-rate). 

Furthermore  a  flat  earth  is  assumed.  In  the  applications  the  use  of  coordinates  different  from  those  may  be 
more  convenient. 


Antenna  Coordinates.  In  ground  coordinates  the  kinematical  state  xp(t)  of  a  platform  moving  at  a  constant 
height  hp  and  a  constant  velocity  rp  =  vp(cos  cpp,  sin  cpp,  0)T  is  given  by  xp(t)  =  (rp(t),rp)T  with  r p(t)  = 
rp(0)  +  rp  t  and  rp(0)  =  (xp,  yp,  hp)r .  Let  the  orientation  of  the  antenna  array  relative  to  the  platform  velocity 
rp  be  defined  by  the  direction  vector  efl  =  (cos  cpa,  sin  cpa,  0)T.  The  cases  cpa  =  0,  n  II  describe  side-  and 
forward-looking  array  antennas,  respectively  (SLAR/FLAR).  Due  to  aircraft  crab,  however,  even  in  case  of 
side-looking  antennas  cpa  ^  0  may  occur  (crab  angle).  The  transforms  xak  =  tg^a[xjf;  tk\,  xjf  =  t a^g[xak,tk] 
between  the  kinematical  states  xsk,  xak  in  ground/antenna  coordinates  at  time  tk,  respectively,  are  given  by: 


t a^g\xsk\tk\  =  Tfl^g  [x8k  -  xp(tk )] 
t  g^a[xak,tk]  =  Tg<_ax£  +xp(tk) 


with  T0<_g  =  Tj_a  and  Tg^a  given  by: 


Lg  <-<* 


Va) 

0 


0 

—  Va) 


COS  (p 

RZ(<P)  =  |  -  sin  (p 
0 


(36) 

(37) 


Sensor  Coordinates.  Let  xak  =  ( r" ,  irak)T  with  rak  =  (xk,  yk,  —hp)T  and  rk  =  (xk,  yk,  0)T  denote  the  kine¬ 
matical  target  state  in  antenna  coordinates  and  let  xk  =  (rk,  cpk,rk,  cpk)T  be  the  corresponding  quantity  in 
sensor  coordinates.  The  non-linear  transforms  between  both  coordinate  systems  are  given  by: 


.  ran  /  aAT7  T~T2  f  yk  xkxk  +  ykyk  Xkyk-xkyk 
t^a[x“]  =  Vxk  +  yk  +  hp,  arctan  — . 


Xk  y/x2k  +yl+h2p  xk  +  yk 


(38) 


ta<-j [x^]  =  Vr2k  -  h2p  ^  cos  cpk,  sin  cpk,  -t 


rkrkco&(pk  .  .  rkrk  sin  cpk  . 

(Pk  sin  cpk,  — - - —  +  cpk  cos  cpk,  0 


r2  -  h2 
rk  nP 


r2  -h2 
rk  nP 


(39) 


The  transfomis  between  ground  and  sensor  coordinates,  C<_g[xf ;  f*],  tg^s  [x/c;  tk  ],  result  from  the  concatena¬ 
tion  of  these  transforms.  The  calculation  of  the  related  jACOBlans,  T5<_g[x^]  =  dts<_g/dxjf  and  Tg^v[Xi  |  = 
dtg<-s/dxk,  is  straightforward. 
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2.2.2  Target  Dynamics 

Let  the  kinematical  state  vector  of  a  ground  moving  target  at  time  tk  be  given  by  its  current  position  r®  = 
(xivxkl’4  jT  and  velocity  r®: 


<  =  (4A)T  =  (x 


fc;l  ’ 


■4;6)T- 


(40) 


Due  to  the  small  agility  of  ground  targets,  acceleration  components  arc  omitted.  In  the  previous  scenario 
(Figure  7)  we  in  particular  have:  x®.3  =  x®.6  =  0. 


Cartesian  Coordinates.  In  the  ground  coordinate  system  the  target  dynamics  is  modeled  by  a  linear  system 
equation  with  additive  white  GAUSSian  noise.  With  the  quantity  Z/;|/(_|  =  vf(l  —  e~2('tk~tk-^/e<)1/2  (discussed 
below),  3D  diagonal  matrices  J  =  diag[l,  1.0],  O  =  diag[0, 0, 0],  and  Vk  ~  N{ 0,  J),  a  realization  close  to 
the  model  in  the  first  talk  (subsection  2.2.1): 

4  =  Fi\k-A-i  +  Gl\k~^k  (4D 

with  matrices  | /c „ ,  and  G®^^  given  by: 

*Vl  =  (i  e-k-irwj).  <42> 

The  modeling  parameters  vt  {limiting  speed)  and  0t  {maneuver  correlation  time)  in  ~Lk\k-\  have  a  clear 
meaning.  In  contrast  to  the  dynamics  model  used  in  [9],  the  model  naturally  introduces  a  ‘speed  limit’  vt, 
while  0t  may  characterize  different  target  types. 


Sensor  Coordinates.  With  the  non-linear  transforms  tg_v,  ts<_g  between  ground  and  sensor  coordinates,  a 
first-order  Taylor  expansion  around  x®_|,fc_]  and  \k-\\k-\  (the  state  estimates  at  t/(-]  in  ground  and  sensor 
coordinates,  respectively,  using  all  associated  measurements  up  to  and  including  tk-i)  yields  a  linearized 
system  equation  in  the  moving  sensor  coordinate  system: 


X/t  =  +  Gk\k-\Vk  +  Uk\k-i 


(43) 


where: 


Ffc|fc-i  -  Ts^g[xAr|A;_1]F/t|/t_1Tg^1.[x<:_i|fc_i],  Gfc|£_i  -  Ts^g[x/t|A;_1]G^|A:_1 

—  1  —  F/c |/c  —  1  [X/c — 1|A:  —  1  [  ,  ^k\k— 1  s[^A:— l|fc— Is  tk—l] 

1  [xfc— i|fc— 1]  —  t^gfF^  ^  _1tg^s[x/t_i|/t-i ,  4-iL  hG- 


2.3  GMTI  Sensor  Model 

In  the  subsequent  considerations  Zk  =  {z‘k  }”10  denotes  a  set  of  /i/{  sensor  reports  (frame  of  observations), 
which  are  detected  at  time  tk  (revisit  time).  In  this  notation  let  z°k  be  the  event  that  at  time  4  the  sensor 
possibly  produced  no  valid  target  detection. 

2.3.1  GMTI  Chacteristics 

Due  to  the  physical  and  technical  reasons  previously  discussed,  the  detection  of  ground  moving  targets  by  air¬ 
borne  radar  is  limited  by  strong  ground  clutter  returns.  This  can  be  much  alleviated  by  STAP  techniques  [13]. 
The  characteristics  of  STAP  processing,  however,  directly  influence  the  GMTI  tracking  performance. 
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Main-Lobe  Clutter.  Let  us  consider  a  ground  patch  located  at  the  position  rf  in  ground  coordinates 
(i.e.  rf  =  o),  which  is  illuminated  by  the  main-lobe  of  the  moving  STAP  radar  sensor.  In  antenna  co¬ 
ordinates  let  r  and  <p  denote  slant  range  and  azimuth  angle  of  the  ground  patch  related  to  the  Cartesian 
position  vector  rac  =  V r2  -  lrp  (cos  cp,  sin  cp,  0)T.  According  to  Equation  36  in  the  moving  antenna  coordi¬ 
nate  system  the  state  vector  of  the  ground  patch  is  given  by  x"  =  (rac,  r")T  with  a  non-zero  velocity  vector 
r"  =  -vp(cos  (pa,  sin  (pa,  0)T.  The  radial  velocity  of  the  ground  patch  therefore  directly  results  from  equa¬ 
tion  39: 


rc(r,  V,  vp,  hp,  cpa)  =  -vp  cos (q>  -  <pa)  V 1  -  ( hp/r )2.  (44) 

Besides  the  location  of  the  ground  patch  described  by  r  and  cp  the  radial  velocity  rc  depends  on  the  parameters 
platform  speed  vp,  platform  height  hp,  and  the  orientation  of  the  antenna  with  respect  to  the  flight  direction 
<Pa- 

For  physical  and  technical  reasons  the  signal  processing  for  airborne  GMTI  radar  [13]  is  unable  to  sep¬ 
arate  a  ground  moving  target  from  the  surrounding  main-lobe  clutter  return  if  the  target’s  radial  velocity  is 
equal  to  rc,  i.e.,  if  its  kinematical  state  vector  x/<  =  {rk,  (pk,  i'k ,  (Pk)T  of  the  target  obeys  the  relation  nc(xk)  =  0 
where  the  function  nc  is  given  by: 


nc(xk)  =  nc(rk,cpk,fk-,vp,hp,cpa )  (45) 

=  fk~  fc(rk,<pk-,vp,  hp,  <pa).  (46) 

In  other  words,  the  equation  nc(x  k)  =  0  defines  the  location  of  the  GMTI-specific  clutter  notch  of  the  sensor  in 
the  state  space  of  a  ground  target  and  as  such  reflects  a  fundamental  physical/technical  fact  without  implying 
any  further  modeling  assumptions. 

Qualitative  Discussion.  An  adequate  modeling  of  this  phenomenon  is  important  to  ground  moving  target 
tracking.  Besides  being  simple  enough  to  be  mathematically  tractable,  the  GMTI  sensor  model  must  reflect 
the  following  qualitative  conclusions : 

1 .  The  detection  probability  of  the  sensor  depends  on  the  kinematical  target  state  and  the  sensor/target 
geometry:  PD  =  PD{xk). 

2.  Pn(xk)  is  small  in  a  certain  region  |«c(x* )|  <  MDV  around  the  clutter  notch  characterized  by  the 
sensor  parameter  MDV. 

3.  Far  from  the  clutter  notch,  the  detection  probability  depends  only  on  the  directivity  pattern  of  the 
sensor  and  the  target  range. 

4.  There  exists  a  comparatively  narrow  transient  region  between  these  two  domains  in  the  state  space  of 
the  ground  target. 

Quantitative  Modeling.  In  a  generic  description  of  the  detection  performance  of  GMTI  sensors  it  seems 
plausible  to  write  Pd  =  P[)(xk)  as  a  product  with  one  factor  reflecting  the  directivity  pattern  and  propagation 
effects  due  to  the  radar  equation  [3],  pd  =  Pn(rk,  (pk ),  the  other  factor  being  related  to  the  clutter  notch.  To 
this  end,  let  us  consider  functions  of  the  following  form: 


"cin-mSk)  1  \ 

1  -  e  2  V  mdv  )  j  (47) 

In  this  expression  the  sensor  parameter  MDV  has  a  clear  and  intuitive  meaning:  In  the  region  defined  by 
|wc(x*)|  <  MDV  we  have  Pd  <  ^  p(].  The  parameter  MDV  is  thus  a  quantitative  measure  of  the  minimum 
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radial  velocity  with  respect  to  the  sensor  platform  which  a  ground  moving  target  must  at  least  have  to  be 
detected  by  the  sensor  (Minimum  Detectable  Velocity).  The  actual  size  of  MDV  depends  on  the  particular 
signal  processor  used. 

For  SWERLING  I  targets  pci  is  given  by:  pci(r,  cp)  =  pF|/[|+sn'(rv)l  plc  false  alarm  probability  pr  and 
the  signal-to-noise  ratio  snr(r,  cp)  =  snro  D(cp)  (o/oq)  (r/ro)~4  according  to  [3].  Let  the  sensor’s  directivity 
pattern  be  described  by  D(cp)  =  sin 2(cp  -  cpa). 

2.3.2  Convoi  Resolution 

Since  in  military  ground  traffic  vehicles  often  move  in  convoys,  at  first  view  resolution  phenomena  seem  to 
be  typical  of  long-range  ground  surveillance.  Due  to  the  asymmetric  effect  of  range  and  angle  resolution, 
however,  Doppler-blindness  in  many  cases  superimposes  resolution  effects:  As  soon  as  convoy  targets  cease 
to  be  resolvable,  they  are  at  the  same  time  buried  in  the  clutter  notch  and  thus  escape  detection.  Vice  versa, 
resolvable  convoy  targets  are  rarely  Doppler-screened.  A  separate  modeling  of  the  sensor  resolution  might 
therefore  be  omitted. 

As  an  example  we  assume  two  targets  moving  in  a  row  along  a  straight  road  with  30  km/h  as  typical  of 
military  applications.  Their  mutual  distance  is  50  m.  The  target/sensor  geometry  is  as  depicted  in  Figure  7. 
Let  the  sensor  resolution  be  given  by:  ar  =  10  m  (range),  n(p  =  0.1°  (azimuth),  ar  =  0.5  m/s  (range -rate). 
Figure  10  shows  the  detection  probabilities  of  both  sensors  (solid  lines).  The  width  of  the  notches  is  larger 
than  in  Figure  8  due  to  the  smaller  convoy  speed.  The  dotted  lines  denote  the  resolution  probabilities  Pr  of 
the  sensors  modeled  according  to  the  discussion  in  the  previous  talk  (subsection  3.2.3): 

pr  —  l  _  e~log2(A r/ar)2  g-log2(A log  2(Ar/«r)2  (4g) 

Ar,  A  cp,  Ar  are  the  distances  between  the  targets  in  sensor  coordinates.  If  P,  is  dominated  by  the  angular 
resolution  (i.e.  Ar  and  Ar  are  small),  Doppler-blindness  occurs.  Outside  of  the  notch  the  high  range/range- 
rate  resolution  guarantees  resolved  returns.  For  an  approach  to  track  ground  moving  convoys  see  [15]. 

2.4  GMTI  Data  Processing 

According  to  Bayes’  formula,  the  processing  of  the  new  sensor  data  Z*  received  at  time  ?/(  makes  use  of  the 
predicted  density  p(xic\Zk~1)  (previously  calculated)  and  the  measurement  likelihood  function  p(Z/- ,  |x*.). 

2.4.1  Likelihood  Function 

The  likelihood  function  piZ^.n^ |x>t )  statistically  describes  what  a  single  frame  of  reports  Z/(  =  {z lk}^L0 
at  time  t *  can  say  about  the  target  state  x* .  Its  particular  structure  is  determines  by  the  current  sensor  data 
frame  and  the  sensor  model  describing  the  sensor’s  properties.  The  likelihood  function  can  be  written  as 
a  sum  over  all  interpretation  hypotheses  regarding  the  origin  of  the  data.  For  well-separated  targets  in  a 
possibly  cluttered  environment  two  classes  of  data  interpretations  can  be  identified: 

1.  The  object  being  considered  was  not  detected,  the  received  sensor  data  are  false. 

2.  The  target  was  detected,  z',  is  the  target  measurement,  all  other  returns  arc  false. 

For  the  sake  of  simplicity  let  us  first  consider  Doppler-unambiguous  measurements.  According  to  the 
previous  modeling  assumptions  we  obtain  up  to  a  factor,  which  is  independent  of  x*  and  the  measurements 
(see  the  discussion  in  the  previous  talk,  subsection  3.3.1),  the  following  expression: 

Hjfc 

p(Zk,nk\xk)<x(l-PD(xk))pF  +  PD(xk)^  J\f  (z£;  Hxfc,  R)  (49) 

n= 1 
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Figure  10:  detection  and  resolution  probability 


with  pf  denoting  the  spatial  false  return  density.  In  this  expression  the  function  nc{Xk)  in  the  definition  of 
the  detection  probability  Poi^k)  (equation  47)  is  linearized  around  the  predicted  state  estimate  x^-i  = 
(rk\k-i,<Pk\k-i,fk\k-i,  (Pk\k-\)T ■  We  obtain: 


Wc(Xfc)  &  Zk~  Hfcxfe  (50) 

with  quantities  Zk  and  H/{,  which  arc  given  by: 


Zk  —  ftc(^k\k— l)  "I"  HfcXfc|£_i 
dnc(xk) 


Ht  =  — 


dxk 


— xfc|fc—  1 


vph2p  cos  cpk\k-\  VpSin^ifc-iV/ 


k\k—\ 


hl 


k\k- 


-1^ 


k\k-\ 


fk\k-l 


(51) 

(52) 

(53) 


The  notation  chosen  indicates  that  the  effect  of  the  GMTI-specific  clutter  notch  on  the  likelihood  function 
can  formally  be  described  by  a  fictitious  measurement  z,k  and  a  corresponding  fictitious  measurement  matrix 
H/(.  The  real  measurements  are  denoted  by  z'k.  By  these  considerations  the  detection  probability  Pn(Xk)  can 
thus  approximated  by  using  a  GAUSSian,  which  linearly  depends  on  the  target's  kinematical  state  vector  x/{ : 


PD(*k)  ~  PD{rk\k-\.<Pk\k-\) 


1 


MDV 
Vlog  2/;r 


^  Hfcx,, 


(54) 


with  a  fictitious  standard  deviation  of  the  measurement  error  related  to  Zk  essentially  given  by  the  minimum 
detectable  velocity  (MDV),  being  a  characteristic  sensor  parameter  of  the  STAP  radar. 
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In  order  to  simplify  the  notation  we  introduce  quantities  z'kJ ,  Ft”7 .  R"7  defined  by: 


i1 = 

H-  = 

r»01  MDV2 

"  2  log  2 

,»0  _  n 

Jk  _  Ac’ 

Hf  =  H, 

R"°  =  R, 

n  —  1, . . 

■,nk 

(55) 

’nk=«-^T- 

H”1  =  (H,Hfc)T, 

R"'=diag[R,^], 

,,nk. 

We  formally  introduce  R00  =  oo.  According  to  the  previous  approximation  the  likelihood  function 
p(Zk ,  nk\xk)  can  therefore  be  written  as  a  weighted  sum  of  GAUSSians: 

nk  i 

p{Zk,  nk\xk)  ct  £  2  *k  ^izk’  R"7)  (56) 

n= 0  y'=0 

with  weighting  factors  given  by  (n  =  1, . . . ,  nk  ): 

=  pf[1  -  pD(rk |fc_i,  cpk |*_i)],  nk  =  pD(rk\k-i,  <pk\k-i)  Pf  ^==  (57) 

xf  =  PD(rk\k-i,  cpk |fc_i),  =  -Pz)(r/t|/t-i,  (P/tp-i)  ;^==  (58) 

2.4.2  Filtering  Process 

According  to  Bayes’  rule  the  processing  of  the  new  sensor  data  Zk  received  at  revisit  time  tk  is  based  on 
the  predicted  density  p(xk\Zk~l)  and  the  likelihood  function  p(Zk,  nk\xk).  With  equation  56  we  obtain  up  to 


a  normalizing  constant: 

p{xk\Zk)  oc  p{Zk,nk\xk)  p(xk\Zk~1)  (59) 

nk  1 

oc  ^  ^  Pk_i  xk  J\f(zk  ;  xk,  R  7)  J\f(xk\  xfe|fc_j,  PJt^_1)-  (60) 

«= 0  i,j= 0 

The  product  of  GAUSSians  in  Equation  60  can  be  rewritten  according  to  the  following  product  formula: 

W(z;  Hx,  R)  W(x;  y,  P)  =  W(x;  Hy,  S)  W(x;  y  +  W(z  -  Hy),  P  -  WSW)  (61) 

where  S  =  HPHT  +  R  W  =  PHTS_1  (62) 


with  compatible  vectors  and  matrices  x,  y,  z,  H,  P,  R.  This  formula  is  a  consequence  of  the  observation  that 
the  left  hand  side  of  Equation  61  can  be  interpreted  as  a  joint  density  p{ z,  x)  =  p(z|x)p(x).  The  right  hand 
side  is  follows  by  computing  p(x|z)  and  p(z).  We  thus  obtain  p{xk\Zk)  as  a  GAUSSian  mixture 

nk  1 

p(*k\zk)  oc  2  2  pT  (x*;  xZ>  FZ)  (63) 

n=0  i,j= 0 

with  mixture  parameters  p„J ,  xkL,  and  P™^  given  by  the  KALMAN-type  following  update  equations: 


=  0 

y1 . 

nij 

-0  % 

(64) 

II 

^3 

1 

O 

O 

nij 

% 

II 

^3 

1 

*■<5. 

•  h"7x'  S'"7  ) 

’  nk  Xk\k-l’  ^k\k-ll 

(65) 

vo/o 

AA:|A: 

—  x( 

_  Ak|k— 1 

nij 

X 

Ak|fc 

-  x'  +  w"'7 

_  Ak|k-1  ^  k\k- 

A<- 

(66) 

p0/0 

17  A;|fc 

-  P' 

_  rfc|fc-l 

p'"/ 

—  p'  _  w"17 

_  rk\k-\  ¥T/c|A:- 

■iS> 

nijT 

k\k-l 

(67) 

■nij 

k\k—\ 

p  /  ij  nj  T  ij 

~  ^klk-l^k  ^k 

■'  o«'7 

Sk 

_  Tjnjni  t T/d"r 

+  Rj, 

-f 

(68) 
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2.4.3  Realization  Aspects 

The  probability  density  p(xic\Zk)  (equation  63)  is  a  mixture  with  4(«/<  +  1)  components.  If  it  is  propagated 
according  to  the  BAYESian  formalism,  an  exponential  growth  of  the  number  of  mixture  components  seems 
to  be  inevitable. 

Following  the  spirit  of  the  techniques  used  in  standard  PDA  or  IMM  methods  previously  discussed,  we 
propose  a  suboptimal  approach  for  keeping  the  number  of  mixture  components  under  control,  which  might 
easily  be  refined  if  required.  Let  us  first  consider  the  following  second-order  approximation: 


where  the  quantities  p™,  p\,  x™,  and  P™  are  given  by  moment  matching. 

The  mixture  components  that  define  the  individual  densities  p,(x^)  are  related  to  the  data  interpretations 
due  to  the  uncertainty  regarding  the  origin  of  the  reports.  They  may  be  handled  by  well-established  mixture 
reduction  techniques  discussed  in  the  previous  talk  (subsection  3.5).  Provided  the  clutter  density  is  low, 
a  second-order  approximation  of  each  p,(x*)  seems  to  be  reasonable  and  results  in  a  GAUSSian  (standard 
PDA).  For  the  sake  of  simplicity  this  approach  has  been  used  in  the  numerical  results  presented  below. 
Hence  the  p(xic\Zk)  is  finally  approximated  by  a  normal  mixture  with  two  components  at  each  time  f&. 
The  generalization  to  more  refined  approximation  techniques  is  straightforward  and  leads  to  MHT  (Multiple 
Hypothesis  Tracking). 

2.5  Discussion 

Figures  11-16  provide  a  qualitative  insight  into  the  effect  of  the  refined  sensor  model  on  target  tracking/data 
fusion.  While  a  high  adaptivity  is  evident  near  the  clutter  notch,  far  from  the  notch  no  difference  to  standard 
filters  is  observed. 

2.5. 1  Effect  of  GMTI-Modeling 

Figures  11a,  lib  display  the  probability  density  functions  resulting  from  processing  the  event  that  a  missing 
detection  occurred  near  the  notch.  To  show  the  most  interesting  features,  the  densities  arc  projected  on  the 
azimuth/  range-rate  plane.  While  the  pdf  of  the  standard  tracker  (Figure  11a)  is  identical  with  the  corre¬ 
sponding  predicted  density,  the  refined  sensor  model  leads  to  a  bimodal  structure  (Figure  lib).  The  broader 
peak  refers  to  the  possible  event  that  the  missing  detection  has  simply  statistical  reasons  as  in  case  of  standard 
filtering,  while  the  sharper  peak  behind  it  reflects  the  hypothesis  that  the  target  was  not  detected  because  it 
is  masked  by  the  clutter  notch.  The  situation  that  the  target  is  buried  in  the  clutter  notch  for  several  revisits 
is  represented  in  Figures  12a,  12b.  Evidently  the  pdf  of  the  standard  filter  totally  faded  away  permitting  no 
reasonable  state  estimation  (Figure  12a).  The  refined  filter,  however,  preserved  a  definite  shape  (Figure  12b). 
This  can  be  explained  as  follows:  Instead  of  actual  sensor  data  the  very  information  that  several  successively 
missing  detections  occurred  was  processed.  This  event  provides  a  hint  to  the  filter  that  the  kinematical  target 
state  probably  obeys  a  certain  relation  determined  by  the  clutter  notch.  Apparently,  this  piece  of  evidence 
proves  to  be  as  valuable  as  a  measurement  of  one  of  the  components  of  the  target  state.  Figures  13a,  13b 
refer  to  the  event  that  a  detection  occurred  near  the  clutter  notch.  While  the  standard  filter  produced  a  simple 
GAUSSian,  the  refined  filter  shows  a  more  complex  structure.  In  fact,  the  pdf  is  a  two-component  mixture 
whose  weighting  factors  differ  in  their  sign  (but  sum  up  to  one).  The  resulting  shape  permits  an  intuitive 
interpretation:  The  sensor  model  inherently  takes  into  account  that  the  target  state  x*  does  not  lead  to  a  small 
value  of  «t(x/{);  otherwise  the  target  would  not  have  been  detected  at  all.  For  this  reason,  the  sharp  cut  in  the 
pdf  simply  indicates  the  location  of  the  clutter  notch. 
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Figure  11:  (a)  standard  filter 


(b)  GMTI  filter  (no  detection) 


Figure  12:  (a)  standard  filter  (b)  GMTI  filter  (within  the  notch) 


Figure  13:  (a)  standard  filter  (b)  GMTI  filter  (detection  near  the  notch) 
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2.5.2  Gain  by  Sensor  Data  Fusion 

Figures  14-16  show  the  probability  densities  of  the  target  position  in  Cartesian  ground  coordinates  after 
filtering.  The  prolated  structure  of  the  pdfs  mirrors  the  predominant  impact  of  cross-range  errors.  Their 
shape  is  rotated  with  respect  to  each  other  due  to  the  different  sensor/target  geometries.  This  effect  can  be 
much  more  pronounced  in  other  situations.  We  indicated  the  true  target  position.  Figures  14a  -  16a  refer  to 
a  regular  tracking  situation  (after  10  min,  see  Figures  7,  8).  Doppler-blindness  occurred  for  sensor  2  during 
the  previous  revisits.  The  probability  densities  shown  in  Figure  14b  -  16b  have  been  calculated  at  a  time 
when  the  target  has  stopped  for  3  min  (Figure  8).  Evidently  in  Figures  14b,  15b  the  dissipation  of  the  density 
functions  is  confined  to  a  particular  direction  according  to  the  GMTI  sensor  model.  Figures  16a,  16b  show 
the  probability  densities  obtained  by  sensor  data  fusion.  In  both  cases  we  observe  a  significant  fusion  gain. 
It  is  a  consequence  of  the  different  orientation  of  the  density  functions  and  leads  to  improved  state  estimates. 
The  result  for  the  stopping  targets  is  particularly  remarkable.  Though  no  sensor  data  are  available  from  both 
sensors,  the  very  fusion  of  the  sensor  output  ‘target  under  track  is  no  longer  detected’  implies  an  improved 
target  localization.  This  is  a  consequence  of  the  different  target/sensor  geometries. 


3  Fusion  with  Context  Information:  Road  Maps 

In  many  practical  cases  even  military  ground  vehicles  move  on  roads,  whose  topographical  coordinates  are 
available  up  to  a  certain  error  (digitized  road  map  information).  In  this  context  it  seems  reasonable  to  describe 
the  kinematical  state  vector  x'k  of  road  targets  at  time  4  by  its  position  on  the  road  4  (i.e.  the  arc  length  of 
the  curve)  and  its  scalar  speed  4:  xrk  =  (4,  4)T.  The  model  for  describing  the  dynamical  behavior  of  road 
targets  is  therefore  a  2D  version  of  equation  41.  By  making  use  of  the  related  transition  density  p(xrk  |x^—1 ) 
the  predicted  density  in  road  coordinates  is  given  by  p(xrk\Zk~l)  =  J  dx[_l  p(xrk |x£_j)  p(x[_  1  \Zk~x)- 

3.1  Modelling  of  Roads 

A  given  road  through  a  real  road  network  is  mathematically  described  by  a  continuous  3D  curve  72*  in 
Cartesian  ground  coordinates.  For  the  sake  of  simplicity  the  effect  of  crossroads  is  not  considered  here.  See 
[21,  22]  for  a  more  detailed  discussion.  Let  72*  be  parameterized  by  the  corresponding  arc  length  /.  The 
exploitation  of  digitized  road  maps  provides  the  data  base  for  a  piecewise  linear  approximation  of  the  road 
curve  TV  :  / 1-»  72* (0  by  a  polygonal  curve  72.  Let  us  furthermore  assume  that  the  curve  72  is  characterized 
by  nr  node  vectors 


Sm  —  72  ( lm ),  tri  —  l,...,  nr. 

From  the  these  quantities  nr—  1  normalized  tangential  vectors 


(70) 


t m  —  (Sm+1  Sm)/l|Sm+l  Sm||,  TH  —  flr  1  (71) 

can  be  derived.  The  EuCLlDian  distance  ||sm+i  -sm||  between  two  adjacent  node  vectors,  however,  is  usually 
not  identical  with  the  distance  Am  =  lm+\  —  lm  actually  covered  by  a  vehicle  when  it  moves  from  sm  to  sm+\ 
along  the  road.  Besides  the  vectors  sm  the  scalar  quantities  Am  >  |sm+i  -  sm|  should  therefore  enter  into  the 
road  model  to  make  it  more  realistic.  The  differences  =  \Am  -  ||sm+i  -  s,„||  |  can  evidently  serve  as  a 
quantitative  measure  of  the  discretization  errors  we  have  to  deal  with.  Using  the  characteristic  functions 

f  1  f*->r  l  G  (On.  6n+l)  nil  /"70\ 

Xm(l)  =  0  dse  .  m  =  0,...,nr,  l0  =  -oo,  l„r+l  =  oo  (72) 

and  s0  =  si,  t0  =  ti,  lo  =  l\,  t„r=t„r_i,  (73) 
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Figure  14:  (a)  tracking 


(b)  target  stop  (sensor  1) 


Figure  15:  (a)  tracking  (b)  target  stop  (sensor  2) 
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we  obtain  a  mathematically  simple  description  of  the  polygon  curve  72,  by  which  the  road  72*  is  approxi¬ 
mated: 


nr 

n  -.i  h-  7 i(i)  =  2  [s»  +  (/ 

m= 0 


with:  72*(/,„)  =  7 2(/m)  =  sm. 


Xni(J) 

m  =  0, . . . ,  nr. 


(74) 

(75) 


3.2  Densities  on  Roads 


The  BAYESian  formalism  previously  discussed  can  directly  be  applied  to  road  targets,  if  it  is  possible  to 
find  a  transformation  operator  7^<_r  by  which  the  predicted  density  p(x'k  |Z/<_I )  in  road  coordinates  can  be 
transformed  into  ground  coordinates: 

Kx*|Zfc-‘)  =  rg^r[p(xrk\Zk~l)\.  (76) 


When  available  in  ground  coordinates,  the  linearized  versions  of  the  transforms  ts<_g  and  tg<_s  (section  2.1) 
can  be  used  to  represent  the  densities  in  sensor  coordinates,  where  the  filtering  step  is  performed.  To  this 
end,  we  write  the  density  p(xKk\Zk~l )  as  a  sum  over  the  nr  +  1  road  segments  considered: 


nr 

p(xgk\Zk~l)  =  2  P(*Sk,rn\Zk~l) 

m= 0 
nr 

=  Yjp(xsk\m,Zk-X)p{m\Zk-1) 

m= 0 

m= 0 

In  equation  79  the  probability 


p(m\Zk~l ) 


dxrkp(m,x[\Zk  1 ) 
dx[  zm(Urx[)  p{x[\Zk~l) 


(77) 

(78) 

(79) 

(80) 

(81) 

(82) 


9-32 


RTO-EN-SET-086 


Tracking  and  Data  Fusion  Applications 


denotes  the  probability  that  the  target  moves  on  the  segment  m  given  the  accumulated  sensor  data  Zk  1 .  The 
matrix  Hr  is  defined  by  Hrx£  =  1^.  Later  on,  it  will  be  intuitively  interpreted  as  a  fictitious  measurement 

matrix.  Since  the  density  p(x'k\Zk~l)  =  Z]=o  p{-i  •^'{xk ;  x'k\k-\  ■  ^k\k-\)  's  a  GAUSSian  mixture  due  to 
the  GMTI  sensor  model,  the  probabilities  can  explicitly  be  expressed  by  error  functions: 

PgU  =  mwjm+1)]  -  mui)]),  m  =  0,...,nr  (83) 


with: 


®(A)  =  1  /  \fln 


rX 


dt  exp(-r/2) 


mj  = 


I  -  H,  x 


rj 

k\k—\ 


Vh,p;;_,h; 


For  the  remaining  term  in  equation  79  standard  probability  reasoning  yields: 

Tg%[p(x'k\Zk-1)\=p(xsk\m,Zk-1) 

dx[  p(xsk,xrk\m,  Zk~l) 

dx[p(xl\xrk,  m)  p(x[\m,  Zk~'). 


(84) 

(85) 

(86) 

(87) 

(88) 


3.2.1  Simple  Roads 

Let  us  first  consider  the  simple  limiting  case  of  a  straight  road  defined  by:  72(/)  =  s  +  It.  Under  GAUSSian 
assumptions  the  transform  from  road  to  ground  coordinates  is  defined  by  the  normal  transition  density: 
t(xLiK+i)  =  ^(4+v  °i)  with  the  affine  transform  tg^r[xr]  =  ( J  (t’)xr  +  (*7/*)  and  am  de¬ 

noting  the  standard  deviation  of  the  mapping  error.  The  transformation  of  the  density  p(xk\Zk~])  into  the 
ground  coordinate  system  is  therefore  described  by  p(xsk\Zk~l)  =  J  dxrk  p(xsk\xrk)  p(x'k\Zk~1)-  The  integra¬ 
tion  can  explicitly  be  carried  out  and  preserves  the  GAUSSian  character  of  the  density  functions  (normal 
mixtures).  The  corresponding  inverse  is  simply  provided  by  a  projection  of  the  density  p(xKk\Zk~l )  on  the 
road.  With  these  transformations  the  previous  considerations  directly  apply. 


3.2.2  Polygonal  Roads 

The  transition  density  p(xsk \xrk,  m)  for  the  road  segment  m  (equation  88)  is  characterized  by  road  map  and  dis¬ 
cretization  errors  (am,  og),  that  may  vary  from  segment  to  segment.  Under  GAUSSian  assumptions  regarding 
the  possible  error  sources,  with  the  affine  transforms  t"Ur[xr]  =  ( *7  t°  )x,.  +  ( )  for  each  individual 
road  segment  m,  and  the  error  standard  deviation  07  =  a~n  +  07,  we  obtain  normal  transition  densities 

p(xg\xr,m)  =  jV(xg;  t"V_r[xr],  of).  (89) 


With  these  preliminaries,  an  application  of  Bayes’  rule  to  the  remaining  density  in  the  integrand  of  equa¬ 
tion  88  yields: 


p(xrk\m,Zk  l)  = 


p(m\xrk)  p(xrk\Zk  !) 

J  dx[  p(m\xrk)  p(xrk\Zk~l) 


(90) 


with  probabilities  p(m\xrk)  given  by: 


p{m\xrk)  =  /m(Hrxrk). 


(91) 
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Up  to  now  the  derivation  was  exact.  Due  to  the  normalization  constant,  however,  the  characteristic  functions 
violate  the  GAUSSian  character  of  the  probability  densities.  To  circumvent  this  problem  we  propose  the 
following  normal  approximation: 


p(m\xrk)  m  exp -  Hrx,.)2/22]  =  km  H rxr,  A2 )  (92) 

with  z™  and  A;n  given  by: 


;2  _ 
Am  - 


—  2 Om+1  T  l m) 
Um+l  —  Ini) 


12 


(93) 

(94) 


The  quantities  zf  and  2  2  can  be  interpreted  as  the  mean  and  variance  of  a  uniform  density  given  by  /„,(/)■ 
From  equation  91  and  the  product  formula  (section  4.3)  we  obtain: 


p(x[\m,Z 


-■>  =  I  pT.,  v(x:.;  X-"/  „  p 


,rmj 


Nfc’  k\k-\ ’  x  k\k 


Li) 


j= 0 


(95) 


pm  j  ffTl  I  O 

with  KALMAN-type  update  equations  for  xfc|fc-1  and  PA.|yt_1,  where  z™,  Afn  aic  in  analogy  to  a  measurement 
and  a  related  measurement  error  variance: 


rmj 

xk\k-i 
p  rmj 
rk\k-l 


=  X 


+  W  mj(z‘ 

k\k—\  + 


m  _  tr 


-  H  xmj  ) 

nrXk\k-\> 


=  VT\L\  -  wrsrw 


(96) 

(97) 


with  “innovation”  covariance  matrices  S,.  and  “KALMANGain”  matrices  Wr  given  by: 

c'nJ  _  rr  p rmj  ttT  ,  ,2 
~  k\k—la'r  '  Am 


r  ^  rk\k-l^r 


Hjsrr 


(98) 

(99) 


The  notation  chosen  indicates  that  the  effect  of  road  map  information  on  the  probability  density  functions 
can  formally  be  described  by  a  fictitious  measurement,  a  corresponding  measurement  matrix,  and  a  fictitious 
measurement  error.  Now  the  integration  in  equation88  can  be  carried  out  explicitly  as  in  the  previously 
discussed  limiting  case.  The  transformation  from  road  to  ground  coordinates  is  thus  known.  In  analogy  to 
the  limiting  case  of  straight  roads,  the  inverse  transform  is  simply  provided  by  individually  projecting  the 
densities  p(xsk\m,  Zk)  on  the  road  (i.e.  after  the  filtering  step  has  been  performed).  Before  the  subsequent 
prediction  is  performed,  it  seem  to  be  reasonable  to  apply  a  second-order  approximation  to  the  mixture 
densities: 


^  p{m\Zk)  p{xrk\m,  Zk) 

m= 0 

(100) 

1=0 

(101) 
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